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The time-dependent stress relaxation for a Rouse model of a crosslinked polymer melt is com- 
pletely determined by the spectrum of eigenvalues of the connectivity matrix. The latter has been 
computed analytically for a mean-field distribution of crosslinks. It shows a Lifshitz tail for small 
eigenvalues and all concentrations below the percolation threshold, giving rise to a stretched expo- 
nential decay of the stress relaxation function in the sol phase. At the critical point the density of 
states is finite for small eigenvalues, resulting in a logarithmic divergence of the viscosity and an al- 
gebraic decay of the stress relaxation function. Numerical diagonalization of the connectivity matrix 
supports the analytical findings and has furthermore been applied to cluster statistics corresponding 
to random bond percolation in two and three dimensions. 



I. INTRODUCTION 

The most striking observation in near critical gels is the anomalous stress relaxation^, which precedes the trans- 
formation of the viscous fluid into an elastic amorphous solid, i.e. the gelation transition. Here, polymer systems are 
considered, where the viscoelastic behavior is controlled by the concentration c of crosslinks connecting monomers 
of different molecules. At a critical concentration c cr i t the gelation transition occurs. Viscoelastic studies by several 
groups have revealed the following characteristic features of stress relaxation: 

• In the sol phase, well below the gelation transition, one observes a stretched exponential decay of the stress 
relaxation function \(t) ~ exp — (t/r*) 13 . 

• The time scale r* ~ e~ z diverges as the critical point is approached. Here e = (c cr ;t — c)/c cr it denotes the 
distance from the critical point. 

• The viscosity 77, which is given as the integral over the stress relaxation function, diverges r\ ~ e~ k as the critical 
point is approached. 

• At the critical point, stress relaxation is algebraic in time: \{t) ~ £~ A . 

Whereas the stretched exponential decay is characteristic of the sol phase and holds for all crosslink concentrations 
c < Ccrit, the last three observations refer to critical behavior, as the gel point is approached. If dynamic scaling holds, 
these findings can be cast in a scaling ansatz for the stress relaxation function \{e.,t), which depends on time and 
distance from the critical point e: 

X {c,t) = e*- k g{t/T*{e)) (1) 

with t* ~ e~ z . Given a certain distance e from the gel point, one expects to see a crossover from an algebraic decay 
at intermediate times to a stretched exponential decay at asymptotically large times. [The scaling ansatz implies 
A = (z — k)/z. Dynamic scaling as implied by Eq. (|l|) is well confirmed experiment allya for the intermediate timg 
regime where xfaf) decays like a power law. However, the values for the exponents scatter considerably. Martin et al.0 
and AdolLet al.El find A = 0.7± 0.05 in agreement with the value A = 0.7± 0.02 of Durand et al.u, whereas Winter and 
coworkers! observe a wide range of exponent values 0.2 < A < 0.9, depending on molecular weight and stoichiometry. 
The experimental support for a universal stretched exponential law is weak. Whereas Martin et al. confirm the 
stretched exponential decay and quote j3 ~ 0.40, other studies reveal non-universal exponents (3. The divergence 
of the time scale r* ~ e~ z in the scaling function was determined in viscoelastic measurements as z = 3.9 ± 0.2cra 
and deduced from static measurements of the shear modulus as z = 4.0 ± 0.6Q. The experimental values for k, the 
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critical exponent of the viscosity, vary in the range 0.7 < k < 1.4. The origin of the scatter in the experimental data 
is not clear. One possible explanation is the size of the critical region, which is known to depend on the degree of 
polymerization. Hence experiments with different samples have to cope with different sizes of the critical region and 
possibly strong crossover effects. 

In this paper we study the simplest dynamic model - Rouse dynamics - in the presence of a time dependent shear 
flow by means of analytical calculations and numerical simulations. Within this model, the frequency dependent 
stress relaxation is completely determined by the spectral properties of the connectivity matrix T, which specifies 
which monomers are crosslinked. As a function of the total concentration of crosslinks c one observes in general a 
percolation transition at a critical concentration c cr it, such that for c < c cr it no macroscopic cluster of connected 
particles exist, whereas for c > c cr it the system percolates. In the context of gelation the fraction of sites in the 
macroscopic cluster has been identified with the gel fraction and the percolation transition has been shown to mark 
the onset of solidificationQ. 

The connectivity matrix T is a positive semidefinite, random matrix, which has been studied in various contexts, e.g. 
diluted ferromagnets, diffusion in sparsely connected spacesa, anomalous relaxation in glassy systems and localizations. 
In all cases the system consists of N vertices (monomers in the context of gelation) which are connected by cN edges 
(crosslinks). In the simplest case (mean field) one chooses the edges independently out of all possible N(N—\)/2 edges. 
The density of eigenvalues can be computed analytically for the above simple distribution and has been discussed in 
Refs. ^,^| in the percolating regime, i.e. c > c cr it- In this paper we focus on the range c < c cr ;t, which corresponds to 
the sol phase and the critical point. For crosslinks of unit strength the spectrum of V is shown to consist of (^-functions 
only. For fluctuating crosslink strength, the spectrum is smooth and shows a Lifshitz tail for small eigenvalues. The 
spectrum determines the time dependent stress relaxation function of Eq. (Q). In mean field theory the exponents 
are found to be = 1/3, A = 1 and z — 3. These results have been confirmed by numerical diagonalisation of the 
connectivity matrix T. 

The latter approach can be extended to finite dimensional connectivities, corresponding to two and three dimensional 
percolation. Qualitatively the results are similar to mean-field theory with, however different exponent values. In 
d — 2 we fpd A « 0.74 and k w 1.19. The latter value is in good agreement with a scalimwelation, which was derived 
previouslyL3 and yields k ~ 1.17 using high precision data for the conductivity exponenttil. Our numcrical-analyis in 
d = 3 is restricted to rather small system sizes (< 20 3 ) as compared to high precision data of Gingold et al£3, because 
these authors restrict their investigations to the conductance, whereas we aim at the density of states, which requires 
a much larger numerical effort. Consequently our exponent values are of relatively low precision as compared to Ref. 
|l2| . Specifically we find A ss 0.83 and k w 0.75, whereas scaling together with high precision data yield for the latter 
A; ~ 0.71. 

The paper is organised as follows: In the following section (Sec. ||) we introduce the dynamic model and the 
observables which we want to discuss and which can be related to the spectrum of eigenvalues of the connectivity 
matrix. In Sec. Ill we present the analytical calculations for the mean-field distribution of crosslinks. We brie 



review the derivation of a selfconsistent equation for the spectrum, which was previously given by Bray and Rodger; 
We then go on to discuss the appearance of Lifshitz tails for small eigenvalues. For crosslinks of unit strength the 
spectrum is shown to consist of a countable set of (5-peaks. We present an analytical scheme to systematically compute 
the spectrum by iteration. We also consider crosslinks of fluctuating stength, for which the spectrum is continuous 
and can be obtained by standard numerical means from the selfconsistent integral equation. In Sec. |y| we present 
results from a numerical diagonalisation of random connectivity matrices. We first compute the spectrum for a mean- 
field distribution of crosslinks and compare it to the analytical results. Next, cluster distributions of random bond 
percolation in 2 and 3 dimensions are considered. Data for the stress relaxation function is presented as well as finite 
size scaling plots for the static shear viscosity. We summarize our results in Sec. [v|. Some detailed calculations have 
been deferred to appendices. \—\\~\ n 

Our paper is an extension of previous work, in which we discussed the static shear viscosityE3£3 and self diffusion^ 
in the sol phase as well as at the gelation transition. There it was shown that the long time limit of the incoherent 
scattering function is determined by the zero eigenvalues of the connectivity matrix, and the static shear viscosity is 
determined by the trace of the Moore-Penrose inverse of the connectivity matrix. Here we focus on the full spectrum 
of eigenvalues, which also determines the decay of the stress relaxation at finite times. 



II. MODEL AND OBSERVABLES 



We consider a system of N identical Brownian particles, each characterized by its time-dependent position vector 
R-i(t) (i — 1, ... ,N) in ci-dimensional space of volume V, i.e. with density p — N/V . M permanent crosslinks are 
introduced between randomly chosen pairs of particles (i e , i' e ), resulting in a crosslink concentration c = M/N. These 
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crosslinks are modelled by a harmonic potential 



, M 

e=l 

whose overall strength is controlled by the parameter a > 0. We use units of energy such that k^T = 1 and allow for 
fluctuations in the strength of crosslinks by introducing the parameter A e . Crosslinks of uniform strength correspond 
to all A e = 1 . In general each crosslink e is assigned independently a random strength A e according to the distribution 
p(A). The connectivity of the particles is specified by the connectivity matrix 

M 

r«/ = ^2 ^ (Siie - S H' e ) (&i'ie - <W e ) ' ( 3 ) 
e=l 

in terms of which the potential reads U — (d/2a 2 ) J2fi'=i ' ^i'- As usua l denotes the Kronecker symbol, 

i.e. Sij = 1 for i = j and zero otherwise. 

We consider purely relaxational dynamics in an externally applied space- and time-dependent velocity field u" xt (r, t): 

dtR?(t) = -i |^(t) + w^R^f),*) + £?(*)• (4) 

Here, Greek indices indicate Cartesian co-ordinates a — x, y, z, . . . , and we will always consider a flow field in the 
^-direction, increasing linearly with y, i.e. 

vZct( r ^) = S o t ,xK(t)r y , (5) 

with a time-dependent shear rate n(t). The noise ^ has zero mean and covariance (^f(t) £,it{t')) = 2C _1 S a ^ 5i^> S(t—t'), 
where 6(t) is the Dirac 5-function. Here, the bracket (• • • ) indicates the average over the realizations of the Gaussian 
noise The relaxation constant is denoted by C- The probability distribution of crosslink configurations Q = 
{i e , igjglias well as the probability distribution of crosslink strengths will be specified later. 

In Ref. M we have computed the shear viscosity in the sol phase for (macro-) molecular units of arbitrary internal 
connectivity. It was shown that the dependence on the crosslink concentration and in particular the critical behavior 
near the gelation threshold is the same for all (macro-) molecular units, as long as we consider identical units with 
a finite degree of polymerization. We expect the same universal behavior for stress relaxation on long time scales, 
which are much larger than the longest internal time scale of a single (macro-) molecule. Hence we specialize to the 
simplest units, namely Brownian particles. 



Relaxation of shear stress 



We aim at the computation of the intrinsic shear stress a a p{i) as a function of the shear rate n(t). For the simple 
shear flow (^|), a linear response relation 

a xy (t) = f drx(* -t)k(t) (6) 

J — oo 

is valid for arbitrary strengths of the shear rate n(t). The linear response or shear relaxation function x(t) is given in 



terms of the connectivity matrix as explained in detail in Ref. 13 

N 



\(./.) = ^>;((l- £„(£)».*!><! -^r.c^) =: £ Tr ((1 - E Q {Q)) exp {-^{Q) \ ) ■ ! < ! 



The matrix Eq denotes the projector onto the subspace of zero eigenvalues of T0. For a time independent shear rate, 
n(t) = Kp-the stress tensor is time independent and related to the shear rate a = prjn via the static shear viscosity, 

^mn^^gT- (8) 



given byE2l 
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Selfdiffusion 



To discuss selfdiffusion we set the externally applied velocity field to zero and focus on the incoherent scattering 
function 



S(q,t) : = T lim ^l^exp{iq(R i (t + T)-R i (T))}^ (9) 
and the squared time delayed displacement 

C(t):= lim (if^R^ + TVR.enA. (10) 



We note that Hi(t + T) — R;(i) is a Gaussian Markov process whose distribution is in the limit T — > oo characterized 
by a vanishing mean and the covariance function 

G w {t) := lim (rRi(t + T) - R(T)] [Rj,(t + T) - R,(T)] ' 



Performing the integral in (Ql]) leads to 



The matrix T is non- negative by inspection (see Eq. (Q)), as it should be to ensure relaxation to equilibrium. The 
scattering function as well as the time delayed displacement can be expressed in terms of Gw (t) via 



1 N 

s(^t) = -J2^p{-q 2 Gu(t)} (13) 



1=1 

and 

N 

X 



1 N 

cw^E G -w- ( 14 ) 



i=l 



Density of eigenvalues 

All dynamic quantities of interest have been expressed in terms of T. Accordingly, once we know the eigenvalues 
{~yi]f = i and eigenvectors of this matrix, we can compute dynamic observables for arbitrary times. In the following, 
we shall discuss the density of eigenvalues 



1 N ~1 

Act (7) = lim - ]T % - 7.) - Jim - Tr S( 7 - T) (15) 

i=l 

for several crosslink distributions. Here ~ denotes the average over crosslink realizations. If one splits off the zero 
eigenvalues, £>tot(7) can be written as 

Aot(7) = ro(c)*( 7 ) + (1 - 7 1 o(c)) J D( 7 ). (16) 

where D(j) is normalized to 1 and contains only the non-zero eigenvalues. If we group the particles into clusters, 
the eigenspace of modes with zero eigenvalues corresponds to vectors that are constant within one clusterEj. In other 
words there is one zero eigenvalue for each cluster and the dimension of the null space is just the number of clusters 
iVd. The weight of zero eigenvalues is simply given by the density of clusters, i.e. T (c) = N c \/N. 
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We restrict ourselves to the density of eigenvalues and do not attempt to compute eigenvectors, which is in general 
more difficult. Hence we can only compute observables which can be written a s j ? J"]^-, (f(T))u, where / is an 
arbitrary function of T. The incoherent scattering function is not of this form (Eq. [13|), whereas the stress relaxation 
function is 

W) = 0-- T ^ c ))pj~ <W7) ex p {-p 7 } • 

The zero eigenvalues are not to be included in the integration, due to the term 1 — .Eo in Eq- (0)- Analogously, the 
averaged viscosity is given by 

, =(1 - ro(c)) ^r d7 32). (18) 

2d Jo 7 

In the same way, the disorder averaged, time delayed displacement is determined by 

d( 7 ) / ^ r 2dt i\ , ml 



C(t) = (l-T {c)) — J d 7 -^ (^1-expj- — 7 |) +T (c)-. (19) 

It can also be expressed as an integral over the time dependent response function 

2 /"* If 

C{t) = ~P Jo dr *( T ) +T °( c V ( 2 °) 



III. MEAN FIELD THEORY 



We consider first the simplest distribution of crosslinks which ignores all correlations between crosslinks, i.e. the 
crosslinks are chosen independently of each other and each pair (i e ,i' e ) of particle indices is realized with equal 
probability. As shown in Ref. [l5| the particle clusters exhibit the analog of a percolation transition at a critical 
crosslink concentration c cr it = 1/2. Below this concentration there is no macroscopic cluster and almost all finite 
clusters are trees. The average number of tree clusters T n with n particles is given in the macroscopic limit by 

lim ^= Tw ^"- 2 (2ce- 2 T (21) 
n^oo N 2cn! y ' 

In particular the total number of clusters per particle is 

T (c) = 1 - c. (22) 

These results are independent of the distribution of crosslink strengths, p(A). 
To compute the density of eigenvalues we introduce the resolvent 



for complex argument Q — 7 + ie, e > 0. In the limit e — > 0, we recover the spectrum from the imaginary part of the 
resolvent according to 

Act (7) = - lim ImG( 7 + ie) . (24) 

7T ej.0 

It can be inferred from Eq. ( |i"5| ) that conversely, -Dtot(7) determines G(Q) via 

G(fl) = /" d 7 ^M. (25) 



5 



A. Disorder average by replicas 

Bray and Rodgersi have shown how to reduce the computation of D to t(7) for crosslinks of unit strength (i.e. all 
A e = 1) to the solution of a nonlinear integral equation. Their derivation is easily generalized to crosslinks of strength 
A which fluctuates according to a given distribution p(X). We restrict ourselves to distributions p(X) such that 



f° dA m 
J —pW <co 



(26) 



holds. It will be shown below (see Eq. (p7[)) that this condition is necessary to ensure a finite viscosity in the sol 
phase. Following Bray and Rodgers we introduce a generating function 



/,.(n^)=p(jp*<^- r «>) 



which determines the resolvent, according to 



G( a) = r im (28) 



The average over the disorder is performed with the replica trick, resulting in 

r ( N n ,, a \ I . N oo N \ 

^= Iin^ ex P U^E^ + ^/ d\ P {\)Ye- lX{ ^ )2/2 -cN . (29) 

J** Vf=io=i^ 27r / V 2 NJo / 

We assume that the connectivity is intensive, limAr^ 00 (c/A r ) = and have introduced the notation <j>i = 
{4>\,4>f, ...,</>") for n-times replicated variables. In the next step one decouples different sites as shown in Ref. || 
and performs a saddle point approximation for large N . This gives rise to a self-consistent equation for a function 
.9°(i) 

^oo / dye [ 

g Q (x) = 2c / dXp(X) J - (30) 

h / d y e (My 2 +2g n (y))/2 



which in turn determines the resolvent according to 



dxx 2 e^ m& / 2+9 (*» 



G(Q) = lim - . (31) 



In the last step of the calculation one assumes a replica symmetric solution for the saddle point equation: 



9 n (£)=9 Q (p) with P=JJ2 x l- ( 32 ) 

The limit n —> can then be performed resulting in the following nonlinear integral equation for g n (p) (cf. Eqs. 
(16,17) in Ref. |) 

g Q (p) = 2e f°dAp(A) exp{-^y 1 

° 1 J (33) 

+ 2ice~ 2c J d\p{\) j dx Xph(iXpx) ex p|~y (P 2 + x 2 ) + l —x 2 + g n (x)\ 
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with <7°(0) = 2c. Here l u (z) are the modified Bessel functions of the first kind. The solution of Eq. (|3^) yields the 
resolvent 



OO J \ ■ />oo 



G(Q) = -J —PW + YcJ dpP9ip) (34) 



and the density of eigenvalues 



Act (7) = ^ limlm dppg^(p)\ 



B. Moments and Lifshitz tails 



(35) 



If all inverse moments M n of the density of non-zero eigenvalues M n := J °° d'-fj n D(j), n € N, exist, one can 
derive the following asymptotic expansion of the resolvent 

G(Sl) = ^ + ^ + c V n» M n+1 (36) 

n— 1 

by expanding the denominator in Eq. ( p5| ) in a geometric series. As we show in Appendix A, the lowest moments are 
given explicitly by 



Mi = 1 
4c 



hi I — !— 1 2c 

1 — 2c 



OO 



dA 2dry 
A (a 



= ^ (37) 



and 

M^-MMm-i^^tlff-l^S (38) 
240c 2 V 1 -2c/ 30c(l-2c) 3 24c(l-2c^ 2 



with P„ := / °° dAA - ™ p(A). Here and in the following we assume that very weak crosslinks are unlikely to occur, 
because we are interested in the small eigenvalues which are due to the geometry of the clusters and not due to the 
appearance of weak links. More precisely we require 

hm — p--- — > -. (39) 
MO |lnA| 2 

The divergence of the moments, Mi and M2, suggests a Lifshitz tail of the density of states of the following form 

D( 7 )c<expj-( 70(1 ~ 2C H \, 7 10, c<|, (40) 

since for positive k, this ansatz implies for the inverse moments 

M„ oc (1 - 2c)- 3 ("- 1 \ cti (41) 

Bray and Rodgers have given a heuristic argument in favor of the ansatz ( fI7j| ) with k = 1/2. They argue that out 
of all clusters for given n the linear one has the smallest eigenvalue, namely 7min = 70 n~~ 2 . There is just one linear 
cluster for given n, so that its contribution to the spectrum is 



D lh 



l£( 2ce <->)\( 7 -£)~e-V^. (42) 



Arbitrary finite clusters may be attached to the chain without altering the dependence of the smallest eigenvalue on 
the length of the chain. If a finite cluster of mass nij is attached to site i of the linear chain, the smallest eigenvalue 
is 7 mm = Replacing 

mj ~m=5^ = -J_ (43) 

E„ T n 1 - 2C 
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leads to jmin — 7 °^I 2c ^ ■ The number of clusters contributing to D(-/) for small 7 is much larger, if attachments are 
taken into account: The probability to find a chain of length n, regardless of attachments, is given by (2c) n . Hence 
the density of eigenvalues is estimated as 

Here, we have expanded ln(2c) ~ 2c — 1 for c sufficiently close to its critical value c cr u = 1/2 to obtain the Lifshitz 
tail near criticality. In Appendix A. 3 we derive rigorous upper and lower bounds for D{^), which prove that D("f) 
has indeed a Lifshitz tail of the form D(j) ~ exp (— y/ h(c) / 7) . We are unable to derive the dependence of h(c) on 
crosslink concentration c, which is however suggested by the lowest order moments (|3^) , (frij) . 

In the following two subsections we shall discuss two special choices for p(X). In the first case all crosslinks are of 
unit strength, giving rise to a point spectrum. In the second case the strength of the crosslinks fluctuates according 
to p(\) = exp (— 1/A)/A 2 . The integral equation ( |33"|) simplifies considerably for this distribution and allows for a 
solution by iteration. 



C. Exact solution of the integral equation for uniform crosslink strengths 

For crosslinks of unit strength, the integral equation (|33|) reduces to Eq. (16) in Ref. |J 

g n (p) = 2ccxp(^ p 2 ) jl + i e-^J^dxphiipx) exp (fcll x * + g ^ j 
= 2cexp(^p 2 ) jl - 2ce~ 2c J™dx J^^exp^O - l)j 2 + 9^) 



(45) 



The second equality follows from a substitution x — > px and from the basic relation between the Bessel functions of 
the first kind J„ and the modified Bessel functions I„, in particular l\{x) = —ii\(ix). 

To get some feeling for the spectrum of eigenvalues, we first consider the case of small c. We then have predominantly 
small clusters, i.e. single particles, dimers, trimers etc. The connectivity matrix of a dimer has eigenvalues \\ = 
0, A2 = 2. A linear chain of 3 particles has eigenvalues {0,1,3}, a linear chain of four particles has eigenvalues 
{0,2,2 + y/2, 2 — V2} and a star with three legs has eigenvalues {0, 1,4}. These are the only trees up to C(c 3 ). 
Hence in this order the spectrum consists of ^-functions at the above eigenvalues, with each cluster contributing to 
the weight of the ^-functions according to its frequency of occurrence. Next we show that ^-functions in the spectrum 
correspond to Gaussian functions g a (p). The ansatz 

g n (p) = 2caexp(-iz(fi)p 2 /2) (46) 

where z — z(Ct) is an arbitrary function of = 7 + ie with Im{2;} < for e > 0, leads to G(Q) = —1 + -, In the limit 
Im{z} — ► each zero ji of Re{z(7i)} = gives rise to a 5 function in the spectrum 

Aot(7) = aE uT7fr\i ( 4? ) 
*f \dz/d^i)\ 

Next, we construct an approximation to the integral equation ( [45| ) by successive iteration. We start with 

So 00 : = 2c. (48) 

The first step of the iteration gives 

5 ?(p) = 2cexp (l-e- 2c r dxJ^e^f'-in-l^+g^-))} (49) 



2 / I Jo W V2 V ' P 2 JU > 

i n 

"2Q-V 
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since the integral on the rhs can be calculated exactly£3. The spectrum consists of a S- function at 7 = 0, 1)1(7) = $(l)- 
The next step of the iteration gives 



.g 2 °(p)=2 C exp(-^ 2 J|l-e- 2c y o dx 3 1 (x) exp I ^(Q - 1)^ + ) ) ( VI i 



2cex P ( -if) |l - e"- Jf ^ J^exp (^(0 - D^j g^exp ) f ^ 



(2c) fe ( i m x 2 



by Taylor expansion of the exponential of ffp(-). Again, the integrals appearing in Eq. ( |52| ) can be computed exactly, 
yielding 



fJ2 



fi (p) = 2 C ^a«exp(-^4V) (53) 



2 



fc=0 

Q (2) ._ £ -2c 2 (2). = fl+ 1 1 (54) 



Note that Ylk=o a fc^ = !• I* 1 this iteration, the spectrum is given by 



1 —2c 00 /o \k 1 

D2M = -^-*(7) + E «- 2C ^b)! 5 ^ " *>" (55) 



Next, we consider a general ansatz for of the form 

L 



9? 



(p) = 2cf>W exp {-^iV}, (56) 



with y^^L aj^ = 1. L is an arbitrary positive integer and will be let tend to 00 below. We insert the ansatz (|5^) into 
Eq. ( ffq ) and use a similar Taylor expansion as above to obtain 

When we now let L — ► 00, we get the expression 

fl&iW - 2c £ agj 1 ' exp (-^V) (58) 



with 



_ -2c TT ^ CCt fc 



fc=0 



and 



4'.*" = 1 + n - 1 - Efao't4" 

We use the notation (/&) to denote a whole sequence of non-negative integers, while Ik (without parentheses) denotes 
the k-th element of the sequence. Out of all possible such sequences we only need those with a finite number of 

W 

non-zero elements in (Ik). The set of all sequences with a finite number of non-zero elements is denoted by {(h)}- 
The summation in Eq. (^8|) thus goes over a countable set and therefore, gf + i(p) is of the same functional form as 

gf'(p)- It is easy to see that J2{(i k )} a (i^P = ^ l 10 ^ a ls° f° r the next iteration. 



non-zero elements. This is because — > as k — > 00, and thus HfeLo " — = if there were infinitely many 
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Since g\ }(p) is an expression of the form of Eq. (|56|), it follows by induction that all g\ (p), i > 2, are of the 
same form. This observation enables us to write down fix-point equations for the coefficients a and the exponential 
prefactors z: 



l (h) 



and 



Ji+l) 



(61) 



As shown in App. [B|, these equations can be solved if the indices on the left and right hand sides are matched by 
mapping the sequence (Ik) that appears as index on the left hand side onto a simple number n = ^2 k lkM k with some 
positive integer M . Afterwards, M is taken to infinity. In the process, a new structure of the coefficients a and z 
emerges: each pair of coefficients (a^, z^) falls into one of infinitely many "classes" of increasing complexity. The first 
three classes are given by the following expressions (the upper index denotes the class), the general form can be found 
in the appendix: 



-2c 



0\n 



1 _ -2c M) 



a (h) = e 



fi ( 2 "4)' fc 



fc=0 



n - i 

Q - 



n - i 



J2k=0 ^ kZ k 



(62) 
(63) 

(64) 



Note that the expression for a higher class automatically contains all of the lower classes as well if the lower-class 
expressions are recursively inserted, e.g. a 



1,0,0,. 



-2c (Zca}) 1 
1! 



a\. 



This remains true in the general case. For higher 
classes, the indices become more complicated, e.g. for class 3 it is necessary to use (1^)) as index on the left hand 
sides. As a shorthand, however, it is convenient to use the notation (Ik) or just k even for the higher classes. It is 
then understood that k itself may stand for a more complicated object like a nested sequence. See App. [b| for details. 
We mention the result that s m , the sum over all a from classes to m, is given by 



„-2c 



'(h 



-2c 



n« 



exp{-2c(l - s™" 1 )} , and 



(65) 
(66) 



As long as c < |, the corresponding fixpoint equation, s — e _2c ( 1 ~ s ) ; has a stable fixpoint at s = 1, which implies 
lim m ^oo s m = 1, as it should be. The quantity 1 — s m is therefore a measure for the quality of an approximation that 
only goes up to class m. We can conclude that for small c only a few classes are sufficient whereas for c close to i 
considerably more are needed. For c > ^, the fixpoint becomes unstable, indicating that the iteration does no longer 
converge to the full solution of the intergral equation due to the appearance of the infinite cluster. 



Implications for the density of states 
Making use of the solution just constructed, the resolvent can be written as 

G(Q) = -1+ lim V-^. 



m— >oo * — ' Zi 

k « 



(67) 



Here, inclusion of a's and z's from classes lower than m in a™ and z™ has been implied as explained above. Analogous 
to Eq. (p7|), this results in the exact density of states 



Act (7)= lim E°™E 



*(7-7S) 

ia*r/07(7K)r 



7 V rfci^ 



that is, a sum of J-peaks located at the roots 7^™ of the respective ^(7) with weight factors 

be proved with Cauchy's integration theorem applied to (z7f >) _1 and Eq. ( |6~I| ) or the more general expression from 



(68) 



It can 



App. pi that J2 i 



Him) 



(h) 

1 holds for every 2™. This property guarantees that the total weight of all peaks 
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in the spectrum is 1 (recall that the sum of all a's is also 1). There is no continuous part of the spectrum, but this 
would change for c > \ due to the appearance of an infinite cluster. 

It is impossible to find the roots of all z m but classes and 1 can be solved exactly. We deduce from Eq. ( |63| ) that 
the roots of z\ are located at 7 n>1 = and 7„.2 = n + 1. The weight factors are easily computed as for the peak 
at and for the peak at n+1. The density of eigenvalues including class and 1 then reads 

^(7) = — SM + E^L*-*). (69) 



Note that this is different from the result of the second iteration, Eq. (p5|), although it contains the same peaks. 

Another consequence of the exact solution of the integral equation is that the density of states does not show scaling 
behavior with respect to c, i.e. it can not be written in the form -Dtot(7) ~ f{lll*{ c )) with some typical 7*(c). This 
follows from the fact that the positions of the peaks are given by the roots of the z's which are independent of c 
and only the weights of the peaks depend on c. This can obviously never result in an exact scaling form: if scaling 
were valid, a small change of 7* would result in a small shift of the peak positions, but they must stay fixed. It will 
be shown below for the fluctuating crosslink strengths that numerical solutions for the eigenvalue density indicate 
that not even an approximate scaling relation holds. This view will furthermore be supported by the results of the 
numerical diagonalization of random matrices for different types of systems. 

To conclude the discussion of the density of states for uniform crosslink strengths, the spectrum from the iterative 
solution of the integral equation is compared with results from numerical diagonalization of T (for details see Sec. IV 
below). Fig. |l| shows the numerically computed spectrum for c = 0.1. Note that there is a peak at 7 = 1, which is not 
present in Eq. (|69|). This "missing peak" can only be found in higher classes, e.g. in £0,1,0,... = 7 2'i 7 5 ~2^g~2 3 2 1 ■ Other 
roots that can easily be identified with peaks in the numerical results are at 2 ± \/2 (stemming from z\ 1 ) or at 

I ± (stemming from Zq 2 )■ Fig- @ shows a direct comparison between the same numerical simulation and a few 
explicitly calculated peaks from classes up to class 3. The agreement regarding the position of the peaks is excellent 
but some weight is still missing from some of the peaks. This weight is expected to be found in higher classes and/or 
in different z's which happen to have a root at the same position. 

D. Numerical integration for special p(X) 
The integral equation (^) simplifies considerably for a special choice of p(X) , namely 

p( A ) = J^cxpj-Ij, (70) 
implying P n = nl. Inserting the ansatz g n (p) =: fn(p 2 /2) into Eq. ( pi| ) leads to the following representation 

i r°° 

G(Q) = -1 + - dxf n (x), (71) 
2c Jo 

where fn(x) is the solution of the ordinary differential equation (see Appendix A. 2 for details) 

fa{x) = -ixf&(x) + 2cexp{-2c + inx + f n (x)}, / n (0) = 2c. (72) 

This allows one to write down the general term in the asymptotic expansion of G(f2) for small f2. Close to the critical 
point the lowest order moments are explicitly given by 

M^ljln^)-^}, c ^i, (73) 

M2 = T5(I^ + 60(^ (74) 
M3= 240(l 4 l2c)a + 105(l-2c) 6 +O « 1 - 2c ^ °A (?5) 
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and 



5762 1159 mtl , 7 , 1 

M ^ 6435 (l-2c) 9 + 72O72O(l-2c) 8 +0((1 - 2c) } ' C ^2' < 76 > 

giving additional support to the conjecture about the Lifshitz tail Eq. (flO|). 

For a numerical evaluation of G(Q) it is more convenient to rewrite Eq. (|35| ) in the following form, 

ff n (p) = 2cV2ipK 1 (V2ip) 

d7?Ki(V2ir/) exp|y?7 2 + 3 (7/)|, 

since in this representation, the integrands do not depend on p and the numerical integration thus needs to be done 
only once per iteration, resulting in time and memory requirements only of the order of the number of integration 
grid points. This allows for high precision computations of g n (p), G(Q), and D(-f). 

Figures || and ^ show the results for the density of eigenvalues from a numerical integration of Eq. ([77]) using a 
Pade approximation in order to extrapolate = 7 + ie to e = 0. There are several noteworthy points to be seen in 
these figures: 

First, we expect to see Lifshitz tails for all c, < c < 1/2, for asymptotically small 7. Precisely at the critical point 
D(~f) goes to a constant as 7 — * 0. For crosslink concentrations close to the critical one, we expect to see a crossover 
between an approximately constant region at intermediate 7 to a Lifshitz tail at very small 7. Since small values of 7 
are hard to access numerically, this crossover makes it difficult to observe the Lifshitz tail, except possibly for small 
c. For intermediate c the data in Fig. || can be described approximately by a straight line but with a slope different 
from — i . This property will be confirmed by the results from the numerical diagonalization presented below. 

A second remarkable point is that the density of states as seen in Fig. || is clearly not suited to a scaling ansatz. 
There are (at least) two different scales contained in the plot: the first is the drop-off length 7 (c) which describes 
the scale on which D(j) goes to for small 7, the other is the position of the maximum, 7 max (c). While 7 goes to 
for c — * i, j max evidently does not; these two features together are obviously incompatible with a scaling ansatz of 
the form f (7) ~ f{l/l*{ c )) with some typical 7*. This finding is in agreement with the observation from the exact 
solution for uniform crosslink strength where scaling was not possible either. Here, however, the statement is even 
stronger than in the previous case since even an approximate scaling relation is ruled out. Note the peculiar feature 
that a second maximum appears in D(p() for small 7 at the percolation threshold c= |. This is not an artifact and is 
confirmed by the numerical diagonalization as shown below. It may even indicate the presence of a third scale since 
the emergence of a maximum can already be suspected in the curves for smaller c. 



E. Stress relaxation 



The characteristic features of the spectrum as discussed above, have important consequences for the stress relaxation 
function. In particular the Lifshitz tail in the spectrum gives rise to an anomalous long time decay of the stress 
relaxation function in the sol phase for all c < c cr j{. The true asymptotic behavior of D(pf) ~ exp(— yh(fi)H), which 
is proven rigorously in Appendix A. 3, implies x(i) ~ exp(— (t/r*)P) with /3 — 1/3. However we are unable to estimate 
the timescale needed to reach the asymptotic regime. For smaller times, the stress relaxation function is characterized 
by effective exponents, just like the spectra in Figs. can be fitted to Lifshitz tails with effective exponents that 
depend on crosslink concentration c. 

The divergence of the timescale T*(e) ~ e~ z is determined by the function h(c). The expansion of the resolvent for 
small fl suggests z — 3. At the critical point the density of eigenvalues is constant as 7 — > 0, implying a logarithmic 
divergence of the static shear viscosity and ~ t~ A with A = 1. 

The absence of scaling in the density of states is also relevant for the stress relaxation function. The presence of 
more than one characteristic scale for the eigenvalues implies more than one characteristic timescale for the stress 
relaxation function. As a consequence, the stress relaxation function does not scale either. This point will be discussed 
further below in the context of numerical diagonalization of the connectivity matrix. Attempts to scale data for the 
time dependent stress relaxation function fail, see Fig. [19. 
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IV. NUMERICAL DIAGONALIZATION 
A. Numerical methods 

In this section the eigenvalue densities D(j) of three different types of random networks are studied numerically: 
mean- field (MF) networks as well as two- and three-dimensional simple square/cubic grids. For the first case, crosslinks 
are allowed for all pairs i, j 'of nodes while for the other networks only crosslinks between neighboring nodes may appear. 
For the finite-dimensional grids we apply periodic boundary conditions in all directions. The size of the networks is 
denoted by N, with N = L d (d = 2,3) for the finite-dimensional cases. For the numerical treatment, we consider 
random graphs with a fixed number M of vertices, i.e. the crosslink concentration is c = M/N. Every crosslink has the 
same probability of occurrence. For the implementation of the graphs on the computer, the LEDA libraryEH was used. 
Network sizes up to N = 10000 (MF), N = 3136 (d = 2) and N = 4096 (d = 3) where studied. For each system size 
up to 10 4 different realizations of the disorder were considered (1000 for the largest sizes). Different concentrations 
of the crosslinks between_Q and the percolation threshold c cr it were treated, where c cr ;t(MF) = 1/2, c a n(d = 2) = 1 
and Cc rit (d = 3) «0.7464t!l 

We consider the same two cases regarding the strength of the crosslinks as above: Either all crosslinks have the same 
strength A = 1 or their strengths are distributed randomly with the probability density givepin (|7Q) . Numerically, 
the random values for the strengths of the crosslinks are drawn using the inversion methocO. A random number 
r is drawn which is uniformly distributed in [0,1]. Then the values of A := — 1/hxr are distributed according to 
(f70|). For testing purposes also some systems were studied where the strengths were uniformly distributed in the 
interval [0.5,1.5]. In all cases no significant deviations of the measurable quantities for different distributions could 
be observed. The main difference is that for crosslinks of unit strength the distribution D(fy) of the eigenvalues is 
dominated by a sum of delta-peaks below the percolation threshold while for crosslinks of continuous strength the 
distribution D(-y) is purely continuous (see below). 

The numerical method works as follows: Random networks are created, witluconstant or random crosslink strengths 
as needed. Then, for each graph the connected components are determinedE3. For each connected component the 
connectivity matrix is calculated, whichis a real symmetric matrix. Therefore, for determining its eigenvalues the QR 
algorithm and the Householder methodEil can be applied. Next, the eigenvalues are sorted in increasing order. Each 
connected component has one smallest eigenvalue 0. Because of numerical errors usually the smallest eigenvalue is 
not zero but quite small, depending on the distribution of the strengths of the crosslinks. Consequently, the smallest 
eigenvalue is assigned the value zero. Finally, the eigenvalues of all components are collected, sorted again, and stored 
for further evaluation for each realization of the network. 



B. Results for the mean- field system 

First, we consider the density D(-y) of non-zero eigenvalues for the mean- field network at the percolation threshold 
c = 1/2. Data for the case p(X) = 8(X — 1) have already been presented in Fig. [j]. Here we consider the case where 
the strengths of the crosslinks are distributed according Eq. (|70|). In Fig. [|the resulting density is shown for different 
system sizes together with the analytical result (obtained from the numerical solution of Eq. ([77])). It can be seen that 
the size N = 10000 is already sufficient to find the analytical behavior for a large range of eigenvalues. Especially the 
"dip" near 7 = 0.15 is validated by the numerical data (see inset). Because of the finite system sizes, arbitrarily small 
eigenvalues cannot be found, thus the numerics disagree from the analytical result in that region. Later we will see 
that this raises several problems when studying dynamical properties of the networks. Nevertheless, the analytical 
result lim 7 _>o-D(7) > can indeed be confirmed by extrapolating the numerical data to the infinite system size. 

The behavior of D{pf) for different crosslink concentrations c is presented in Fig. ^. Once more, the numerical 
(N = 10000) and the analytical results agree very well. For small 7, the logarithm of the density should behave as 
~ — 7 -1 / 2 (Lifshitz tail). Fig. [7] shows the logarithm of .0(7) in a double logarithmic plot in complete analogy to the 
solution of the integral equation (Fig. |^). Presumably, the system size of N = 10000 is still too small in order to 
observe the asymptotic behavior of the density of states for small eigenvalues. 

C. Results for the finite-dimensional systems 

Now we consider three-dimensional systems, which are believed to describe real polymer networks more appropri- 
ately. The density of eigenvalues for the case where all crosslinks have the same strength, p(\) = S(X — 1), is shown 
in Fig. p| for TV = 16 3 and c = 0.2. Similar to the mean-field case, a collection of delta-peaks is obtained. Since this 
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kind of distribution is more difficult to analyze, we turn again to the model where the strengths of the bonds have the 
broad distribution ((7^). Results for the largest system size N = 16 3 and different crosslink concentrations are shown 
in Fig. 0. Below the percolation transition c cr it ~ 0.7464 the distribution exhibits a maximum and converges to 
for small eigenvalues, similar to the mean- field case. But at the transition D("f) diverges when the zero eigenvalue is 
approached (see also inset). Below we will show that this behavior changes the divergence of the viscosity near the 
percolation threshold. The eigenvalue densities for the two-dimensional network look qualitatively similar and are 
therefore not shown here. The true asymptotic behavior as 7 — > is difficult to access, just as in the mean-field case. 

The changes in the spectrum as compared to the mean-field case also effect the stress relaxation, that we investigate 
next. First, the viscosity given by 

r? = (l-T (c)) f°^d7 (78) 
Jo 7 

is considered. Here, irrelevant prefactors have been omitted for simplicity, cf. Eq. (|l^) for the complete expression. 
In the numerical calculation wc compute 

7.>0 ' 

for each realization and subsequently average over different realizations of the disorder to obtain rj. Whereas for the 

mean- field network the viscosity diverges logarithmically for c — ► c cr it, for finite-dimensional systems a divergence 

77(c) ~ ( c cr it — c)~ fc is expected. The reason for the different divergences is the manner in which D(-f) behaves for 

small 7 at the percolation threshold: for the mean-field network, lim 7 ^o -^(7) is finite, but for the finite-dimensional 

grids D(^) diverges as 7 — ► 0. The exponent k, which describes the viscosity near the percolation threshold, can be 

hi 



determined, similar to the usual finite-size scaling relational for the percolation transition, from 

V (c,L)=L- k ^((c-c clit )L x ^), (80) 

where rj is a universal function and v is the exponent describing the divergence of the correlation length when 
approaching the percolation transition. The use of finite-size scaling enables us to circumvent the problems which are 
posed by the lack of very small eigenvalues of finite graphs. 

By plotting r\L k l v against (c — c C rit)L 1 < v with correct parameters v and k the datapoints for different system sizes 
for c ps Ccrit should collapse onto a single curve. We have taken the values v (d = 2) = 4/3 and v (d = 3) = 0.88 from 
the literatures and adjusted k/v. The best collapse near c cr it was obtained with k(d = 2) = 1.19 and k(d = 3) = 0.75. 
The results are presented in Figs. [l(] (d — 3) and |ll| (d — 2). The values we have obtained for the different distributions 
of the crosslink strengths agree within the error bars. 1— . 

The value of k for two dimensions agrees very well-with the result k ~ 1.17 found previously by Broderix et al.EJ, 
using the high precision simulations-oLGingold et alB3. The result for the three-dimensional case (k — 0.75)is slightly 
worse in comparison with k ~ 0.7lE3li3. The reason is that here only small system sizes up to 20 3 could be treated 
due to the fact that all eigenvalues are calculated. If one is only interested in k, it is computationally less expensive to 
compute the Moore-Penrose inverse of the connectivity matrix. Thereby one might be able to study system sizes as 
large as those used in Ref.li3. For the realizations treated here, we have checked other characteristic results concerning 
the percolation transition, like the critical exponent cr, which describes the behavior of the cluster-size distribution. 
The finite-size scaling plots have a poor quality for this quantity, too, resulting in a rather low precision of the exponent 
values. Additionally, we have observed a systematic drift in our results: By including even smaller system sizes, the 
scaling plot results in k = 0.89 which differs even more from the value obtained before. Consequently, we believe 
that larger system sizes are needed, to obtain a more reliable result for k via numerical diagonalization of random 
connectivity matrices. 

Next, the behavior of the stress relaxation function (again omitting irrelevant prefactors and using dimensionless 
thne$-t) 

X (t) = (1 - T (c)) / £>( 7 ) exp(- 7 i) d 1 (81) 



is investigated, see Eq. (|17j ) for the complete expression. The functions were obtained by first calculating D(j) and 
then numerically integrating it. It would take too much time on the computer to first calculate x(i) for each realization 
by directly summing up the contributions and then average over the disorder. Here, we have investigated the systems 
with the continuously distributed crosslink strengths because they result in continuous eigenvalue densities where it 
is easier to obtain stable numerical data. 
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In Fig. |12| the numerical results for the mean-field network, the d = 2 and the d = 3 model for the largest sizes 
(c = c cr it) are shown. As mentioned before, the numerical simulations are restricted to finite sizes of the networks 
and to a finite number of realizations of the disorder. Therefore, the eigenvalue densities .0(7) always have a smallest 
eigenvalue 7 m i n with D(j) — for 7 < 7 m i n . Consequently, the long-time behavior is dominated by an exponential 
decrease, exp(— 7mini)> irrespective of the true form of x(t)- This results in a negative curvature in the double- 
logarithmic plot for long times. Thus, in the numerical results, the asymptotic form of the relaxation function is 
visible only for intermediate times (see Fig. |l2|). At c = c cr it a x(~t) ~ t A behavior is expected. By fitting we obtain 
A = 1.029(5) (mean field), A = 0.830(2) (d = 3), and A = 0.741(2) [d = 2). The results for the mean-field case 
is known exactly to be A = 1. The discrepancy comes again from the finite sizes of the networks: Indeed, we have 
observed that for smaller networks a value of the exponent is obtained which is even larger. So the result A = 1 seems 
to be confirmed. The value for the three-dimensional grid is compatible with the large range of results obtained in 
experimental. 

The behavior of x(t) f° r different concentrations c of the crosslinks is shown in Figs. [D| (mean-field) and |bl| (d = 3). 
In both cases we find exponential decay for the longest times due to finite system size. For intermediate times a 
stretched exponential behavior x(i) ~ exp(— (t/r) 13 ) is visible. At least for finite system sizes the exponent j3 seems 
to be non-universal, we find values ranging from = 0.5 for small crosslink concentrations down to (3 = 0.2 close to 
the percolation threshold. We suspect that the accessible times are too short to see the true asymptotic behavior, 
which at least in mean-filed theory is known to be a stretched exponential with exponent (3 = 1/3, resulting from the 
Lifshitz tail in the density of states. For small times %(i) decreases like t~ A and x(0) = 1 by definition. 

Moreover, this variation of the exponent (3 makes it impossible to observe a scaling form %(i) ~ t~ A g(t/r), were r 
is a typical time scale which diverges like r ~ (c cr ;t — c)~ z when approaching the percolation threshold. For the mean- 
field network, the expectations from the Lifshitz tails are z — 3 and A = 1, while g(t) is the stretched exponential 
function, but it was already mentioned that there seems to be no scaling possible due to the existence of more than 
one scale. In Fig. [l5| a scaling plot of is shown. x(t)t A is plotted against t x (c cr ;t — c) 2 for mean-field networks 
having different concentrations c. We have used only the regions below the finite-size asymptotic behavior (j3 = 1). 
It can be seen that the quality of the collapse is rather bad, explained by the variation of (3 with c. One might think 
that near the transition c ~ c cr it the scaling may be better. But there the collapse is even worse (not shown), because 
even larger systems are necessary to reach the asymptotic regime for the small eigenvalues, as explained before. 

For the finite-dimensional systems, the quality of the scaling-plot is similar. Therefore, it is not possible to make a 
reliable estimate for the dynamical exponent z in that case. 

V. CONCLUSIONS 

Within our model, the dynamics of a crosslinked polymer melt is determined completely by the eigenvalue and 
eigenvector spectrum of the connectivity matrix T. In this paper we have focused on some properties which are 
determined by the eigenvalues alone (e.g. the stress relaxation function) since the eigenvectors are hard to obtain. We 
have used three different methods to examine the eigenvalue spectrum: First, the construction of an exact solution 
for the averaged eigenvalue density for a fixed crosslink strength, second, a very precise numerical solution for the 
case of varying crosslink strengths, and third, a numerical diagonalization of random connectivity matrices. 

The first method allowed for some exact results regarding the eigenvalue spectrum. It could be shown that the 
eigenvalue spectrum consists of a very complicated but countable set of 5-peaks, some of which could be calculated 
and compared with results from numerical diagonalization. Furthermore it showed that the eigenvalue density does 
not show (exact) scaling behavior. 

The second model of fluctuating crosslink strengths has the advantage that the eigenvalue spectrum becomes a 
continuous function instead of an inscrutable sum of <5-peaks. Additionally, it allowed for a fast numerical integration 
scheme. From these numerical solutions it could be inferred that the expected Lifshitz-tail behavior for small 7 only 
seems to set in for extremely small 7, smaller than were accessible numerically. For this reason, the stress relaxation 
function does not show a stretched exponential form with exponent j3 = i within the accessible time window. 
Instead, for the times that could be reached, there seems to be a regime where an apparent stretched exponential with 
a crosslink concentration-dependent and thus non-universal (3 can be observed. Furthermore, in numerical evaluations 
of the eigenvalue spectrum again scaling could not be observed, not even approximately, since at least two, possibly 
three or more, different 7-scales with different c-dependence could be identified. As a consequence, the stress relaxation 
function does not scale cither. 

The third method, numerical diagonalization, confirmed all results obtained so far very well. In particular it showed 
that the stress relaxation shows stretched-exponential behavior with a concentration-dependent exponent (3 and it 
showed the failure of scaling of the stress relaxation function. It confirmed, however, the experimental findings that 



15 



at the critical concentration, the stress relaxation function decays algebraically with exponent A. For the mean-field 
model, both theory and numerics yield the exponent A = 1. Furthermore, numerical diagonalization allows for going 
beyond the mean-field approach. Results were obtained for connectivity matrices on two- and three-dimensional cubic 
lattices. Unlike the mean-field case, the density of eigenvalues now diverges at the critical concentration as 7 — * 0. 
Apart from this difference, the other results are qualitatively very similar to the mean-field case except that the 
exponents A and (3 have different values. 
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APPENDIX A: LOW FREQUENCY EXPANSION OF THE RESOLVENT 



1. General p(A) 

The low frequency expansion is derived from an alternative form of the integral equation (|33|). We start from Eq. 
( ^p| ) and recall an integral representation of the n-dimensional Laplacian (see also Eqs. (3.47-51)) in Ref. |l3| 



'2 

(Al) 

p=|x| 



We use this representation in the numerator of Eq. ( |30| ) and take the limit n — > 0. To evaluate the denominator of 
Eq. ( ^0| ) we observe that 

lim fdxfn(\x\) =/o(O) + 0(n). (A2) 

J 

Both steps taken together lead to the following self-consistent equation for g n (p) 

1 / d 2 d M f m 



P°(p) = 2c e -y o dAp( A)exp|-(^-— jjexpj^WWj, (A3) 

which is of course equivalent to the integral equation (|33|), but much better suited for a low frequency expansion. 

To that end we rescale variables according to x — s/Qp and ty a {x) — g n {x/\fTi). The self consistent equation then 
reads 

*"<*> - 2 "" 21 f dA "« «-{^ (a? " is)} a A? 2 + " (A4) 

We look for a solution in terms of a power series in VL 

* n (x)=Y / m^ ] (x). (A5) 

3=0 

The resolvent can then be expressed in terms of &j (x) as follows 

G(Q) = -Pi + — ^ W dxx^>, 



(x) (A6) 



with P n = J °° dAA " p(X) as defined after Eq. (|38|). The lowest order term obeys the equation 

*o(a;) = 2cexp (-2c) exp ( ^x 2 + ^ a (x) J , (A7) 
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which is solved by 

*o(z) = -W (-2cexp (-2c) exp {ix 2 /2)) . (A8) 
Here W denotes the principal branch of Lambert's W^-function, defined as the solution of 

W(x) exp (W(x)) = x. (A9) 



From Eq. ( A7) one derives the following property of the lowest order solution 



which allows for an exact computation of the integral 



1 f°° d 

dxx^oix) = -- dx— (l-^ (x)f = 2c(c-l). (All) 

2 ,/n dx 



The next two terms are given by 



1 A /^ 2 

1 - #o(x) 2i \dx 2 xdx 



^)= , , T , — ^(3-2 ~ — l^oW (A12) 



T , , -1 / d 2 d \ fP 2 P? \ ( d 2 d . 

* 2{x) = T^(x) {d^-^){Y + T^kd ) x~dx) * o{x) - { AU 1 



The integrals J dxx^j(x) can be performed similar to Eq. (All), using the properties of Lambert's W^-function. The 
computations, however, become increasingly tedious, so that higher order terms have only been computed for the 
special distribution p(X), see below. 



2. Special p(X) 



We start from Eq. (A3) and introduce the abbreviation D p := -jy — yj^,- For the special choice p(X) = jz exp (— 1/A), 
one can perform the average over p(X) analytically 



dX 
X 2 
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(l+iDg/2) 
A 



exp 
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iD, 
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in 



(A14) 



The resulting differential equation simplifies, if we introduce the function fn(p 2 /2) := g (p) 

fn( P 2 /2)+ip 2 /2f£(p 2 /2) = 2 C exp(-2c)exp( l n P 2 /2 + /^( j0 2 /2)). (A15) 

Introducing the new variable x = p 2 /2 leads to the differential Eq. ([72]) quoted in the main part of the paper. For the 
low frequency expansion it is convenient to introduce yet another variable, y = tlx, in terms of which the differential 
equation for hn{Clx) :— ffi{x)reads 



hn(y) ~ iy^h'^y) = 2c exp (-2c) exp (iy + h n (y)). 



(A16) 



The ansatz hn(y) — YlTLaity'' hj(y) then yields 



h n (y) - iyK-i(y) = h (y) 



l d n 



■exp [j^inyhM 



(A17) 



The left hand side is linear in h n (y), so that Eq. ( A17 ) is easily iterated. 
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3. Proof of the existence of a Lifshitz tail in D(i) 



The aim of this appendix is to prove that the density of eigenvalues D(j) shows a Lifshitz-tail behavior for 7 — > 
and c < 1/2. For the proof, it is convenient to make use of the eigenvalue distribution function, F(j) := jZ^ D(j'). 
This can be done without loss of generality because if ^(7) has a Lifshitz tail, so does D(-f). It will be shown that 
F(-f) lies between two bounds which, taken together, assert the Lifshitz behavior. 

For a given realization of a system with N vertices (or polymers), the corresponding -F/v(7) can be written, using 
a decomposition into the K clusters of the realization, 



k T , .. K 



" ' N N 

k=l k=l 



where is the connectivity matrix of the fc-th cluster and Eq is the projector on the nullspace of Tk- In the 
macroscopic limit, N — > 00, this yields 

00 

F( 7 ) = (1 - c)6( 7 ) + J2 T " ( Tr [(! - E o )e(7 - r(T„))]> (A19) 

n=l 

due to self-averaging. The bracket (• • • ) means averaging over the set of all numbered trees {T n } of size n of which 
there are n n ~ 2 . T(T n ) denotes the connectivity matrix corresponding to the tree T n . The average number of trees of 



size n per vertex is denoted by r n and is given byliS 

n n— 2 1 
= " ( 2ce - 2c )" = -J_ n - 5 / 2 e-" h ( c )-/(")/" (A20) 
2cn! V ; 2cv / 2¥ V 7 

according to Stirling's formula with h(c) = 2c — 1 — ln(2c) and some function f(n) with < f(n) < 1. 

The smallest non-zero eigenvalue of r(7^,) is certainly greater than or equal to the smallest non-zero eigenvalue of 
the linear cluster with n vertices, which is proportional to n~ 2 , i.e. 

r(T n ) > 4 (A21) 

(except for the zero eigenvalue) with some a independent of n. This results in 

Tr[(l - ££)6( 7 - r(T n ))] < (n - 1)6( 7 - c0i 2 ) (A22) 



F( 7 ) < 1 - c+ ( n _ 1 ) T »' for T > °' ( A23 ) 

n>\/ a/7 

For 7 — > 0, the sum can be approximated by an integral, 

F(7) < 1 - c + = / n- 3 / 2 e-"' l(c) (A24) 



2cv27T Jy/a/f 

"'-'^q'"-!-""^! <A25 » 

This is the lower bound for f (7L 

For the upper bound, Eq. (A19) will be used again. Explicitly, one has for 7 > 

1 00 1 

F( 7 ) = 1 - c + - -^- 2c T E Ml - W(7 - r(T w ))] (A26) 

"=1 ' {T„} 
1 °° 1 

> 1 - c + - £ ^(2ce" 2c )" £ Tr[(l - ^)6( 7 - r(£\))] (A27) 

C n=i n - {£ X } 
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where the inner sum has been restricted on the set of linear numbered trees There are n!/2 such linear trees, 

such that 

oo 

F( 7 ) = 1 - c + ± e-" (,l(c)+1) Tr[(l - £ ")6( 7 - r(£\))]- (A28) 

n— 2 

Next, the trace, which is a sum of non-negative terms, is estimated by just one of the terms. In particular, Tr[(l — 
£?o)0( 7 — > 9(7 — ot/\ e ), corresponding to the smallest eigenvalue of C\. This finally gives 

F( 7 )>l-c+i J2 e~ n{h{c)+1] (A29) 

n~>-\J a /f 

« 1 - c + exp |- (1 + .(c))} (A30) 

for the lower bound. 

The upper and the lower bound together imply 

hm HM^-l + c)! = i (A31) 
7^0 I ln-y| 2 v ; 



or, even stronger, 



Vah(c) < - lim 7 1/2 ln(F( 7 ) - 1 + c) < y/a(h(c) + 1), (A32) 
7— >o 



which is the sought-for Lifshitz tail behavior. 



APPENDIX B: DETAILS OF THE EXACT SOLUTION OF THE INTEGRAL EQUATION 

1. Solution of the integral equation 

It is not obvious how to solve the fix-point equations (|6l|), because the coefficients of the i-th iteration are labeled 
by an index, and a subsequent iteration gives rise to coefficients which are labelled by a sequence (Ik)- We therefore 
try to map the sequence (Ik) that appears as index onto a number by writing n := X)fcLn IfiM with some MeN. 
For this to be a one-to-one map, we need to restrict all Ik to be < M. This restriction will be removed later when we 
let M — ► 00. The sequence (Ik) can be reconstructed from n by writing n in the number system of base M. Let this 
be indicated by Ik — (n)^f ■ 

The fix-point equations can now be written down as: 

k=0 { >k ■ 
1 

The equations for a n can be solved independently from those for z n . We start with a n . Successively solving the 



z n = l + 7^ \ V^oo 7TW7 ■ ( B2 ) 



system of equations (Bl) by inspection gives 



a = er' c (B3) 
. 2c (2cooT for 1 < n < M (B4) 



a„ = e 



n 

M-l 



TT ^ff: for M < n < M M (B5) 
to 

M-l M-l / \(«)f f i„fi 

n n ( fc ° +Mfci+ (n ;r' lfcM r 1 forM^<n<M^ M ( B6 ) 

fc =0 fe M -i=0 \ )ko+Mk\—' 
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and so on. The coefficient ao is obviously independent of all other a n . This property will be called "class 0". 
a\, . . . , q,m-i only depend on ao: this will be termed "class 1". Analogously clm, ■ ■ ■ , &m m -i are m class 2 as they 
only depend on a's from classes and 1. 

Now we can let M tend to infinity. Classes and 1 are simple (the upper index now denotes the class) : 



2c (2ce- 2c r 



n > 1 



(B7) 
(B8) 



If we drop the constraint n > 1, Eq. (BS) automatically contains class 



For the higher classes, as we are now considering M — > oo, indexing via a number n is no longer possible. Instead, 
for class 2, we have to revert to using a finite sequence as index. For class 3, even this is not sufficient and a nested 
sequence (liki)) is needed: 



4 1 - e~ 2c 



n 

k=0 

n 



(2co|) 



length of (l k ) > 1, 



{(fe.)} 



kki) ] 



(B9) 
(BIO) 



If the constraint (length of (lk) > 1) is dropped and if the explicit expressions for the a from the lower classes are 
recursively inserted, all classes up to class 2 are contained in one formula, Eq. (B9). An analogous statement holds 
for Eq. pIcD . 

In general, for class m, the index will be of the form (lr k \) with m nesting levels. The general result is thus 



'(>■;) 



2ca 



(fc 



<fc ) 



" J )) 



n 

)} 



(fc ) 
'('•i) 



(Bll) 



With the same reasoning as above we can calculate the z n . We find the same classes, and the results are: 



n 



tt - 1 



n — viz® 

0,-1- nz^ 



O _ 7 r l 



(B12) 
(B13) 
(B14) 



s (I t * )) 

■fa) 



n-E 



{{* )}H fc . ) 

"'•f«l '>«) 



fi-i-E 



{(k ) z i\ ) 



„m— 1 



(B15) 



2. Properties of the solution 

If one asks for the total weight of a particular peak at, say, some 70 (up to class m), one has to find all finite 
solutions (lk) of the diophantic equation 

00 

7o-E^r~ 1 (^) = °- (Bw) 

fe=0 
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This is possible in some special cases, e.g. for 70 = 1 in class 2. Since z*(l) = 1, it follows that Eq. (B16) is satisfied 
if and only if exactly one entry of (//.) equals 1 whereas all the others are 0. Adding up all of the weights yields 

e _2c ( (2ce _2c — l) e 2ce 2 " + 1 ) as the total weight of S(j — 1) from class 2. 



The have several noteworthy properties most of which are easy to prove by induction over m and are therefore 
listed below without proof. 

1. is a rational function of 7 with integer coefficients. 

2. The degree of the numerator is the same as that of the denominator. 

3. The coefficient of the highest power is 1 in both numerator and denominator. 

4. z^k^ is a strictly monotonically decreasing function (except at its poles). 

5. All roots and poles of are located on the non- negative real axis. 

6. zff\ has exactly as many poles as roots. Roots and poles alternate, starting with a root at 0. 

7. There are exactly one more roots of z™^ than there are poles in X)feLo 

g z m -1 

8. The sum y\ — 57^"(T(z l fe )i) over all roots 7(? fe )j °f z (i k ) ec i ua ls 1- As stated in the text, this can be proved 
using Cauchy's integration theorem. 

Consider now some z™^ and choose (Ik) such that only the n-th entry is nonzero. Then we have 

_ ; m-1 

<-o.«,.o....= 7 7 1 im ;v^ (Bl7) 

Between two of its poles (see the list of properties above), z™ -1 is a continuous function that maps bijectively onto 
the real numbers; therefore there exists a 7™ in this interval such that zff^d™) — 0- Moreover, when l n — > 00, the 

7™ converge to the root 7™" 1 of z™ -1 in this interval. Since z™ -1 is monotonously decreasing, 7™ < 7™~ 1 . This 
implies that for every peak in the spectrum there are infinitely many other peaks to the left of it in any arbitrarily 
small interval around this peak. This also applies recursively for each of these satellite peaks! Only the peak at is 
different: as stated in the list above, all roots of are > and thus there are no satellite peaks of (5(7). 
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FIG. 1. Numerical simulation of the density of states for c = 0.1. 
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FIG. 2. Comparison between the simulation (solid lines) and some selected peaks calculated from the exact solution (dashed 
lines) for c = 0.1. The analytical peaks have been slightly shifted to the right for better comparison, otherwise they would be 
indistinguishable from the numerical peaks. 
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FIG. 3. Density D(^) of non-zero eigenvalues as a function of 7 for p(\) given by (^) for several concentrations c. For 
c = 0, the density is given by D{y) = ^p(^). 




FIG. 4. Double logarithm of the density D(j) of non-zero eigenvalues as a function of In 7 for several concentrations c. 
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FIG. 5. Density 0(7) of non-zero eigenvalues for the mean-field network at the percolation threshold c = 0.5 from numerical 
diagonalization. The solid line is the analytical result, which is hardly distinguishable from the result for N = 10000. The inset 
magnifies the region 7 € [0, 0.4], where the numerical results for the largest system size N = 10000 are shown by circles. 




FIG. 6. Density D(-y) of non-zero eigenvalues for the mean-field network for different concentrations c. The lines are the 
analytical results while the numerical data is shown by the symbols. 
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FIG. 8. Density D(j) of non-zero eigenvalues for the cubic network with all bonds having the same strength A = 1 
(c = 0.2, N = 16 3 ). Similar to the case of the mean-field network, a sum of delta-peaks with strongly varying heights is 
obtained. 
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FIG. 10. Finite-size scaling plot of the viscosity rj(c, L) for the three-dimensional grid. A scaling behavior of 
r](c,L) = L~ k / V rj{{c — c cr it)I/ 1/ '' / ) is assumed. Using v = 0.88 and k = 0.75 the points for L = 10,13,16,20 collapse onto 
one curve near the critical concentration. 
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FIG. 11. Finite-size scaling plot of the viscosity r)(c,L) for the two-dimensional grid. A scaling behavior of 
7]{c,L) = L- k/ "f)((c- Ccrit)!/ 17 ") is assumed. Using v = 4/3 and k = 1.19 the points for L = 10,14,20,28,40,56 collapse 
onto one curve near the critical concentration. Since finite systems are treated, the maximum of rj(c) is below the critical 
concentration c cr i t = 1 of the infinite lattice. 
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FIG. 12. Stress relaxation function x(f) at the the critical concentration c = c cr i t for the three types of models considered 
here, with continuously distributed strengths of the crosslinks in all three cases. Shown are the results for the largest sizes which 
could be treated with sufficient accuracy. For the part of the long-time behavior which is accessible to the numerical simulations, 
a x{t) ~ t" A behavior is visible. From fitting we obtain A = 1.029 (mean field), A = 0.830(2) (d=3) and A = 0.741(2) (d=2). 
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FIG. 13. Rescaled stress relaxation lunction — ln(x(t)t A ) as a function of the time for the mean-field network (A = 1.029) 
with different concentrations c of the crosslinks. The straight lines correspond to stretched exponentials with exponent j3 = 0.332 
and (3 = 1. 




FIG. 14. Rescaled stress relaxation function — ln(x(i)i A ) as a function of the time for the three-dimensional network 
(A = 0.830) with different concentrations c of the crosslinks. The straight lines correspond to stretched exponentials with 
exponent /3 = 0.386 and (5 = 1. 
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FIG. 15. Scaling plot for the stress relaxation xW* A as a function of i(c cr it — c) z for the mean-field network (N = 10 4 , 
randomly distributed strengths of crosslinks) with the values A = 1, z — 3. 
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